LXM: Better Splittable Pseudorandom Number Generators (and Almost as Fast)
En el 2021 Guy Steele y Sebastiano Vigna presentaron un generador de números pseudo aleatorios (PRNG) que se construye a partir de otro PRNG en donde colaboró Guy Steele llamado SplitMix.
La principal propiedad de ambos PRNGs es que estos generadores son partibles (splittable), es decir, se pueden dividir en dos nuevos generadores estadísticamente independientes, lo cual es de gran utilidad en ambientes concurrentes.
LXM se basa en combinar tres ideas: un generador lineal congruente (la L), un generador basado en XORs (la X) y el resultado de la combinación de esos dos generadores utilizarlo como input en una función mezcladora (la M).
En este trabajo nos centraremos específicamente en la parte práctica del algoritmo y lo re-implementaremos en nuestro propio código, evitando hablar tanto de la concurrencia como de la propiedad partible del algoritmo (lo cual consiste de la mayoría del paper de Guy Steele y Sebastiano Vigna).
Las ideas principales de este trabajo surgen de la sección 2 del paper (THE LXM GENERATION ALGORITHM)
def LXM(s, m=2891336453, a=1310709051, x0=11565345, x1=5242836, k=64):
# https://gist.github.com/trietptm/5cd60ed6add5adad6a34098ce255949a
rotl = lambda val, r_bits, max_bits: \
(val << r_bits%max_bits) & (2**max_bits-1) | \
((val & (2**max_bits-1)) >> (max_bits-(r_bits%max_bits)))
while True:
# Combining operation
z = s + x0
# Mixing function (lea64)
z = (z ^ (z >> 32)) * 0xdaba0b6eb09322e3
z = (z ^ (z >> 32)) * 0xdaba0b6eb09322e3
z = (z ^ (z >> 32))
z &= (2**k) - 1
# Update the LCG subgenerator
s = (m * s + a)
# Update the XBG subgenerator (xoroshiro128v1_0)
q0 = x0; q1 = x1
q1 ^= q0
q0 = rotl(q0, 24, k)
q0 = q0 ^ q1 ^ (q1 << 16)
q1 = rotl(q1, 37, k)
x0 = q0; x1 = q1
# Return result (modified to be in interval 0-1)
yield z/2**k
n = 5
print(f"{n} números generados al azar, en el intervalo 0-1:")
g = LXM(42)
for r, i in zip(g,range(n)):
print(f"* {r}")
5 números generados al azar, en el intervalo 0-1: * 0.793781427453092 * 0.2155520740899848 * 0.7433470653358116 * 0.28660629734277526 * 0.515362051632894
scatterplot.show()
En el gráfico de dispersión se pretende ver si la generación de los números sigue algún patrón. En el gráfico realizado no se observa ninguno, ya que no se ven formas claras, como por ejemplo una línea recta, ni zonas con una densidad marcadamente superior a la del resto. Lo que se observa es una "nube de puntos", que es lo que debería verse en una distribución uniforme.
scatter3D.show()
Text(0.5, 0, 'X(n)+2')
histogramas.show()
En un histograma se muestra la cantidad de observaciones de la variable estudiada según el intervalo. Permite observar como se distribuyen los números generados. Por ejemplo, si las barras graficadas tuvieran diferencias significativas en sus alturas esto querría decir que la distribución no es uniforme. En los histogramas hechos, podemos ver que a medida que aumenta la cantidad de muestras usada, y consecuentemente se reduce el tamaño del intervalo, la diferencia que se observa entre las barras disminuye. Esto culmina en el último gráfico, que usa 52016 muestras e intervalos de tamaño 0.06, en el que la longitud de las barras es prácticamente igual. Por esta razón podemos decir que utilizando una cantidad de muestras suficientes (ej: 52016), las muestras tenderán a distribuirse uniformemente en los intervalos.
boxplot.show()
Un boxplot muestra los cuartiles de los números generados, el soporte de la distribución y los outliers. Para el caso actual, se puede utilizar para comparar esos datos con los de una distribución uniforme. Así, se puede ver que:
qqplot.show()
En el QQplot se comparan los cuantiles de la variable a estudiar, en este caso los de los números generados con LXM, con los de una variable ideal, en este caso una distribución uniforme. Este gráfico permite observar si la probabilidad acumulada a izquierda de cada valor de la variable a estudiar se aleja o no del de la variable ideal. Se puede ver que si bien en el primer gráfico se observan claras diferencias entre los cuantiles teóricos y los de la muestra, estas pueden ser atribuidas a la baja cantidad de muestras tomadas para dicho gráfico. Lo anterior se debe a que dichas diferencias desaparecen al aumentar la cantidad de muestras en los gráficos siguientes. De hecho, se puede apreciar que en el gráfico hecho utilizando 1000 muestras, aproximadamente entre 0.1 y 0.4 los cuantiles de la muestra se encuentran apenas por debajo de los teóricos, pero esto se corrige en el gráfico hecho 10000 muestras. Por lo tanto, podemos decir que tomando una cantidad de muestras lo suficientemente grande (ej:10000), los cuantiles de los números generados y los de una distribución aleatoria no presentan diferencias.
Con este método comparamos dos distribuciones de probabilidad. Por una parte, usaremos el generador que implementamos, y por otra parte, para la distribución esperada.
Con este método comparamos dos distribuciones de probabilidad. Por una parte, la observada, que es la que obtenemos utilizando nuestro generador LXM. Por otra parte, la distribución esperada, que es la distribución que queremos corroborar que la observada sigue.
from scipy.stats import chi2
import numpy as np
import math
class Chi2:
def __init__(self, numero_de_muestras=1000, g = None):
self._generador = g or LXM(42)
self._numero_de_muestras = numero_de_muestras
def realizar_prueba(self):
frecuencias_observadas = self.frecuencias_observadas()
f_esperada = 1/10 * sum(frecuencias_observadas)
d2 = sum([(f_observada - f_esperada)**2 for f_observada in frecuencias_observadas]) / f_esperada
limite_superior = chi2.ppf(0.95, df=5)
print(f'Frecuencias: {frecuencias_observadas}')
print('Estadístico: {:.2f}'.format(d2))
print(f'Límite superior: {round(limite_superior, 2)}')
self.validar_hipotesis(d2, limite_superior)
def calcular_d2(self):
numerador = (f_observada - f_esperada) ** 2
denominador = f_esperada
d_al_cuadrado = numerador / denominados # acá es sumatoria
def validar_hipotesis(self, d2, limite_superior):
if d2 <= limite_superior:
print("El test acepta la hipotesis nula")
else:
print("El test rechaza la hipótesis nula")
def frecuencias_observadas(self):
# Estas son las que tengo que obtener
numeros_generados = []
for i in range(self._numero_de_muestras):
numeros_generados.append(next(self._generador))
frecuencias = [0] * 10
for numero in numeros_generados:
pos = math.floor(numero*10)
if pos < 10:
frecuencias[pos] += 1
return frecuencias
def pseudo_chi2():
# Obtengo freuencias observadas de 0 a 9 (enteros)
# La frecuencia esperada es 1/10 por la cantidad de números aleatorios que genero
# Con estos valores ahora entonces calculo d2
# d2 = sum (frecuencia_observada - frecuencia_esperada) ** 2 / frecuencia_esperada
# Finalmente para chi2.ppf(0.95, df=5) vemos si d2 es mayor o menos al limite superior
return
# Corremos nuestro experimento
Chi2(numero_de_muestras=10000).realizar_prueba()
Frecuencias: [949, 998, 1019, 1007, 1008, 998, 1029, 968, 1001, 1023] Estadístico: 5.48 Límite superior: 11.07 El test acepta la hipotesis nula
Este test consiste en contar la cantidad de números aleatorios generados que no pertenecen a un intervalo dado, tal que consigamos 'gaps'. Para este test decidimos generar 1000 gaps y utilizar el intervalo [0.2, 0,5].
Finalmente aplicamos chi2 para validar (o rechazar) la hipótesis.
import matplotlib.pyplot as plt
class Gap:
def __init__(self, cota_inferior, cota_superior, g = None):
self._generador = g or LXM(42)
self._cota_inferior = cota_inferior
self._cota_superior = cota_superior
self._p_0 = cota_superior - cota_inferior
self._cantidad_de_gaps = 1000
def realizar_prueba(self):
# Espero a la primera vez que tengo valor dentro de la cota
gaps = self.generar_gaps()
self.plotear_gaps(gaps)
frecuencias_observadas = self.generar_frecuencias(gaps)
d2 = 0
# Aplicamos chi2
for i, frecuencia_observada in enumerate(frecuencias_observadas):
f_esperada = self.frecuencia_esperada(i)
numerador = (frecuencia_observada - f_esperada) ** 2
denominador = f_esperada
d2 += numerador / denominador
# Validamos hipótesis
k = len(frecuencias_observadas)
limite_superior = chi2.ppf(0.95, df=k-1)
print('Estadístico: {:.2f}'.format(d2))
print(f'Límite superior: {round(limite_superior, 2)}')
self.validar_hipotesis(d2, limite_superior)
def validar_hipotesis(self, d2, limite_superior):
if d2 <= limite_superior:
print("El test acepta la hipotesis nula")
else:
print("El test rechaza la hipótesis nula")
def generar_gaps(self):
gaps = []
# Generamos 1000 gaps
while(len(gaps) < self._cantidad_de_gaps):
while (True):
n = next(self._generador)
if n < self._cota_superior and n > self._cota_inferior:
break
gap = 0
while(True):
n = next(self._generador)
if n < self._cota_superior and n > self._cota_inferior:
gaps.append(gap)
break
else:
gap +=1
return gaps
def plotear_gaps(self, gaps):
plt.hist(gaps, color='b', bins="sturges", alpha=0.8, density=True, rwidth=0.8, align='left')
plt.show()
def generar_frecuencias(self, gaps):
frecuencias = [0] * (max(gaps) + 1)
for gap in gaps:
frecuencias[gap] += 1
return frecuencias
def frecuencia_esperada(self, n):
return (((1 - self._p_0)**n) * self._p_0) * self._cantidad_de_gaps
def pseudo_gap():
# Obtengo gaps entre un rango de valores
# (dentro de [0,1] ya que son los valores que genera el LXM programado
# Con estos valores ahora entonces calculo d2
# d2 = sum (frecuencia_observada - frecuencia_esperada) ** 2 / frecuencia_esperada
# La frecuencia esperada podrá ser distinta para cada valor, ya que
# la misma sigue una distribución geométrica
# Finalmente para chi2.ppf(0.95, df=cantidad_de_gaps) vemos si d2 es mayor o menos al limite superior
return
# Corremos nuestro experimento
Gap(cota_inferior=0.2, cota_superior=0.5).realizar_prueba()
Estadístico: 21.92 Límite superior: 31.41 El test acepta la hipotesis nula
Con este métod vamos a validar que las variables sean independientes.
from scipy.stats import chi2_contingency
class TestIndependencia:
def __init__(self, g = None):
self._generador = g or LXM(42)
self._cantidad_de_vectores = 500
def realizar_prueba(self):
nros_x = []
nros_y = []
for i in range(self._cantidad_de_vectores):
nros_x.append(next(self._generador))
nros_y.append(next(self._generador))
self.plot_vectores(nros_x, nros_y)
d2, p, dof, ex = chi2_contingency([nros_x, nros_y], correction=True)
# Validamos hipótesis
limite_superior = chi2.ppf(0.95, df=dof)
print('Estadístico: {:.2f}'.format(d2))
print(f'Límite superior: {round(limite_superior, 2)}')
self.validar_hipotesis(d2, limite_superior)
def validar_hipotesis(self, d2, limite_superior):
if d2 <= limite_superior:
print("El test acepta la hipotesis nula")
else:
print("El test rechaza la hipótesis nula")
def plot_vectores(self, nros_x, nros_y):
plt.title(f"Dispersión de {self._cantidad_de_vectores} números generados al azar")
plt.ylabel("X(n)-1")
plt.xlabel("X(n)")
plt.scatter(nros_x,nros_y,50)
plt.show()
# Corremos nuestro experimento
TestIndependencia().realizar_prueba()
Estadístico: 87.60 Límite superior: 552.07 El test acepta la hipotesis nula
Utilizando el generador implementado en el ejercicio 1:
gen = LXM(42)
def generar_exponencial(lambda_):
u_1 = next(gen)
return -np.log(u_1)/lambda_
def cociente(t):
return e**((-(t-1)**2)/2)
def transformar_gaussiana(x, mean, std):
return x*std + mean
def generar_normal_0_1():
no_aceptada = True
while(no_aceptada):
exp_1 = generar_exponencial(1)
prob = cociente(exp_1)
if(next(gen) < prob):
no_aceptada = False
z = exp_1
if(next(gen) < 0.5):
z = -z
return z
m1, m2, m3 = [], [], [] #muestras
for i in range(100):
z1 = generar_normal_0_1()
z2 = transformar_gaussiana(z1, 10, 2)
m1.append(z2)
for i in range(1000):
z1 = generar_normal_0_1()
z2 = transformar_gaussiana(z1, 10, 2)
m2.append(z2)
for i in range(10000):
z1 = generar_normal_0_1()
z2 = transformar_gaussiana(z1, 10, 2)
m3.append(z2)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 5))
fig.tight_layout(pad=3)
ax1.hist(m1, 50)
ax1.set_xlabel('Valores generados', fontsize=15)
ax1.set_title("Histograma (n = 100)", fontsize=20)
ax2.hist(m2, 50)
ax2.set_xlabel('Valores generados', fontsize=15)
ax2.set_title("Histograma (n = 1000)", fontsize=20)
ax3.hist(m3, 50)
ax3.set_xlabel('Valores generados', fontsize=15)
ax3.set_title("Histograma (n = 10000)", fontsize=20)
plt.show();
El test de Kolmogorov-Smirnov se utiliza para decidir si una muestra proviene de una población con una distribución específica.
Se tiene una distribución de probabilidad acumulada, la cual se propone como distribución teórica de la muestra, y otra distribución de probabilidad empirica, generada a partir de la muestra.
En este caso tenemos las hipótesis:
m1_ks, m2_ks, m3_ks = [], [], [] #muestras
for i in range(100):
z1 = generar_normal_0_1()
z2 = transformar_gaussiana(z1, 10, 2)
m1_ks.append(z2)
for i in range(1000):
z1 = generar_normal_0_1()
z2 = transformar_gaussiana(z1, 10, 2)
m2_ks.append(z2)
for i in range(5000):
z1 = generar_normal_0_1()
z2 = transformar_gaussiana(z1, 10, 2)
m3_ks.append(z2)
state = np.random.default_rng()
dist_teorica_1 = sp.stats.norm(loc=10, scale=2).rvs(size=100, random_state=state)
dist_teorica_2 = sp.stats.norm(loc=10, scale=2).rvs(size=1000, random_state=state)
dist_teorica_3 = sp.stats.norm(loc=10, scale=2).rvs(size=5000, random_state=state)
figsts, (ax1sts, ax2sts, ax3sts) = plt.subplots(1, 3, figsize=(10, 3))
figsts.tight_layout(pad=3)
yval1 = np.linspace(0.0, 1.0,num=100)
ax1sts.step(x=np.sort(m1_ks), y=yval1)
ax1sts.step(x=np.sort(dist_teorica_1), y=yval1)
ax1sts.set_xlabel('Valores generados (n=100)', fontsize=8)
ax1sts.set_title("Probabilidad Acumulada", fontsize=10)
yval2 = np.linspace(0.0, 1.0,num=1000)
ax2sts.step(x=np.sort(m2_ks), y=yval2)
ax2sts.step(x=np.sort(dist_teorica_2), y=yval2)
ax2sts.set_xlabel('Valores generados (n=1000)', fontsize=8)
ax2sts.set_title("Probabilidad Acumulada", fontsize=10)
yval3 = np.linspace(0.0, 1.0,num=5000)
ax3sts.step(x=np.sort(m3_ks), y=yval3)
ax3sts.step(x=np.sort(dist_teorica_3), y=yval3)
ax3sts.set_xlabel('Valores generados (n=5000)', fontsize=8)
ax3sts.set_title("Probabilidad Acumulada", fontsize=10)
Text(0.5, 1.0, 'Probabilidad Acumulada')
figks, (ax1ks, ax2ks, ax3ks) = plt.subplots(1, 3, figsize=(20, 5))
figks.tight_layout(pad=3)
ax1ks.hist(m1_ks, 10)
ax1ks.set_xlabel('Valores generados', fontsize=15)
ax1ks.set_title("Histograma (n = 100)", fontsize=20)
ax2ks.hist(m2_ks, 25)
ax2ks.set_xlabel('Valores generados', fontsize=15)
ax2ks.set_title("Histograma (n = 1000)", fontsize=20)
ax3ks.hist(m3_ks, 50)
ax3ks.set_xlabel('Valores generados', fontsize=15)
ax3ks.set_title("Histograma (n = 5000)", fontsize=20)
plt.show();
from scipy import stats
ks1 = stats.kstest(m1_ks, dist_teorica_1, alternative='two-sided')
ks2 = stats.kstest(m2_ks, dist_teorica_2, alternative='two-sided')
ks3 = stats.kstest(m3_ks, dist_teorica_3, alternative='two-sided')
display ("Valor del estadístico del test Kolmogorov-Smirnov con n = 100: " + str(ks1.statistic))
display ("P-valor del test Kolmogorov-Smirnov con n = 100: " + str(ks1.pvalue))
display ("Valor del estadístico del test Kolmogorov-Smirnov con n = 1000: " + str(ks2.statistic))
display ("P-valor del test Kolmogorov-Smirnov con n = 1000: " + str(ks2.pvalue))
display ("Valor del estadístico del test Kolmogorov-Smirnov con n = 5000: " + str(ks3.statistic))
display ("P-valor del test Kolmogorov-Smirnov con n = 5000: " + str(ks3.pvalue))
'Valor del estadístico del test Kolmogorov-Smirnov con n = 100: 0.08'
'P-valor del test Kolmogorov-Smirnov con n = 100: 0.9084105017744525'
'Valor del estadístico del test Kolmogorov-Smirnov con n = 1000: 0.061'
'P-valor del test Kolmogorov-Smirnov con n = 1000: 0.04839715079181246'
'Valor del estadístico del test Kolmogorov-Smirnov con n = 5000: 0.0144'
'P-valor del test Kolmogorov-Smirnov con n = 5000: 0.6777877521935483'
Dado que en todos los casos el p-valor es superior a 0.05, no rechazamos la hipótesis de que los datos provienen de una distribución normal.
El test de Shapiro-Wilk está diseñado para decidir si una muestra aleatoria proviene de una distribución normal. En este caso se cuenta con las siguientes hipótesis:
$H_0$: Los datos provienen de una distribución normal
$H_1$: Los datos no provienen de una distribución normal
El test provee como resultados un estadístico y un p-valor. Cuando el p-valor obtenido es superior a cierto umbral --usualmente 0.05--, decimos que no podemos rechazar la hipótesis de que los datos provienen de una distribución normal.
La implementación utilizada (scipy) indica que para muestras de más de 5000 números, se pierde precisión en los valores del test. Por lo que se realizará el test para muestras de tamaños N = 100, 1000 y 5000. (Ver: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.shapiro.html?highlight=shapiro)
m1_sw, m2_sw, m3_sw = [], [], [] #muestras
for i in range(100):
z1 = generar_normal_0_1()
z2 = transformar_gaussiana(z1, 10, 2)
m1_sw.append(z2)
for i in range(1000):
z1 = generar_normal_0_1()
z2 = transformar_gaussiana(z1, 10, 2)
m2_sw.append(z2)
for i in range(5000):
z1 = generar_normal_0_1()
z2 = transformar_gaussiana(z1, 10, 2)
m3_sw.append(z2)
figsw, (ax1sw, ax2sw, ax3sw) = plt.subplots(1, 3, figsize=(20, 5))
figsw.tight_layout(pad=3)
ax1sw.hist(m1_sw, 10)
ax1sw.set_xlabel('Valores generados', fontsize=15)
ax1sw.set_title("Histograma (n = 100)", fontsize=20)
ax2sw.hist(m2_sw, 25)
ax2sw.set_xlabel('Valores generados', fontsize=15)
ax2sw.set_title("Histograma (n = 1000)", fontsize=20)
ax3sw.hist(m3_sw, 50)
ax3sw.set_xlabel('Valores generados', fontsize=15)
ax3sw.set_title("Histograma (n = 5000)", fontsize=20)
plt.show();
from scipy import stats
shapiro1 = stats.shapiro(m1_sw)
shapiro2 = stats.shapiro(m2_sw)
shapiro3 = stats.shapiro(m3_sw)
display ("Valor del estadístico del test Shapiro-Wilk con n = 100: " + str(shapiro1.statistic))
display ("P-valor del test Shapiro-Wilk con n = 100: " + str(shapiro1.pvalue))
display ("Valor del estadístico del test Shapiro-Wilk con n = 1000: " + str(shapiro2.statistic))
display ("P-valor del test Shapiro-Wilk con n = 1000: " + str(shapiro2.pvalue))
display ("Valor del estadístico del test Shapiro-Wilk con n = 5000: " + str(shapiro3.statistic))
display ("P-valor del test Shapiro-Wilk con n = 5000: " + str(shapiro3.pvalue))
'Valor del estadístico del test Shapiro-Wilk con n = 100: 0.9904636740684509'
'P-valor del test Shapiro-Wilk con n = 100: 0.7020654678344727'
'Valor del estadístico del test Shapiro-Wilk con n = 1000: 0.9989786744117737'
'P-valor del test Shapiro-Wilk con n = 1000: 0.8619867563247681'
'Valor del estadístico del test Shapiro-Wilk con n = 5000: 0.9995430111885071'
'P-valor del test Shapiro-Wilk con n = 5000: 0.28750666975975037'
Dado que el p-valor es superior a 0.05, no rechazamos la hipótesis de que los datos provienen de una distribución normal.
import os
wd = os.getcwd()
if 'TP1' not in wd: wd += '/TP1'
arribos_real = []
with open(f"{wd}/tiempos_entre_arribos.txt") as f:
for l in f.readlines(): arribos_real.append(float(l.strip()))
tasa = len(arribos_real)/sum(arribos_real)
print(f"Tengo {len(arribos_real)} arribos en {sum(arribos_real)} horas")
print(f"Mi tasa de arribos empírica es de {tasa} arribos por hora")
Tengo 10000 arribos en 1011.1994386169154 horas Mi tasa de arribos empírica es de 9.889245996494681 arribos por hora
# Para generar eventos de Poisson:
# 1. Generar muestras de una distribución exponencial Z,
# con parametro Lambda: z1, z2, z3...
# 2. Definir tiempos de eventos conteo de Poisson:
# t0 = 0, t1 = t0 + z1, t2 = t1 + z2....
def iter_exponencial(lmbda, generador, n):
"""Genera una secuencia que siga una distribución exponencial,
dado un generador que sigue una distribución uniforme"""
u = np.array([next(generador) for i in range(int(n))])
x = -np.log(1-u)/lmbda
return iter(x)
def generar_poisson(limite, lmbda, generador):
expo = iter_exponencial(lmbda, generador, lmbda*limite*1.5)
arribos = []
t = 0
while t <= limite:
arribo = next(expo)
arribos.append(arribo)
t += arribo
return arribos
limite = 24 * 31 # Un mes
arribos = generar_poisson(limite, tasa, LXM(42))
print(f"Simulé {len(arribos)} arribos en {sum(arribos)} horas")
print(f"Mi tasa de arribos simulada es de {len(arribos)/sum(arribos)} arribos por hora")
Simulé 7300 arribos en 744.064522714941 horas Mi tasa de arribos simulada es de 9.810977109033201 arribos por hora
def comparar_poissones(poissones_con_label, n = 15):
for arribos, label in poissones_con_label:
plt.step(np.cumsum(arribos[:n]),range(n),where= 'post',label=label)
plt.legend()
plt.title(f"Comparación de los primeros {n} arribos")
plt.show()
comparar_poissones([(arribos_real, "Poisson Empirica"), (arribos, "Poisson Simulada")])
limite = 96
poissones = []
for i in range(1000):
poissones.append(generar_poisson(limite, tasa, LXM(i)))
comparar_poissones([(p, f"Poisson Generada #{i}") for i, p in enumerate(poissones[:5])])
import math
# PMF de poisson, a mano, para comparar con los valores teóricos
pmf = lambda x, lmbda: ((lmbda ** x) * (math.e ** -lmbda)) / (math.factorial(x))
print("Probabilidad que el primer vehículo arribe antes de los 10 minutos.")
diez_min = 1/6
# La probabilidad de que haya al menos un arribo en 10 minutos equivale a
# 1 - la probabilidad de que no haya ningun arribo en 10 minutos
print(f"* Probabilidad Teórica => {1 - pmf(0, tasa*diez_min)}")
favorable = []
for p in poissones:
favorable.append(p[0] <= diez_min)
print(f"* Probabilidad Simulada => {sum(favorable)/len(favorable)}")
Probabilidad que el primer vehículo arribe antes de los 10 minutos. * Probabilidad Teórica => 0.8076055651533732 * Probabilidad Simulada => 0.818
print("Probabilidad que el undécimo vehículo arribe después de los 60 minutos.")
sesenta_min = 1
# La probabilidad de que el vehiculo 11 llegué despues de 60 minutos equivale a
# la sumatoria de que hayan llegado 0,1,2,...,10 vehiculos en 60 minutos
# Ojo con los off-by-one! `for i in range(11)` === de 0 a 10
print(f"* Probabilidad Teórica => {sum([pmf(i, tasa*sesenta_min) for i in range(11)])}")
favorable = []
for p in poissones:
arribos_acumulados = np.cumsum(p)
# El arribo en el indice [10] === El undecimo arribo
favorable.append(arribos_acumulados[10] >= 1)
print(f"* Probabilidad Simulada => {sum(favorable)/len(favorable)}")
Probabilidad que el undécimo vehículo arribe después de los 60 minutos. * Probabilidad Teórica => 0.596893339441658 * Probabilidad Simulada => 0.579
# Como nuestra pmf a mano da overflow al calcular un número grande como lambda a la 750,
# importamos directamente el modulo de scipy
from scipy.stats import poisson
print("Probabilidad que arriben al menos 750 vehículos antes de las 72 horas.")
setentidos_horas = 72
# La probabilidad de que arriben al menos 750 vehiculos en 72 horas equivale a
# 1 - la sumatoria de que lleguen 0,1,...,749 en 72 horas
print(f"* Probabilidad Teórica => {1 - sum([poisson.pmf(i, tasa*setentidos_horas) for i in range(750)])}")
favorable = []
for p in poissones:
arribos_acumulados = np.cumsum(p)
favorable.append(arribos_acumulados[749] <= 72)
print(f"* Probabilidad Simulada => {sum(favorable)/len(favorable)}")
Probabilidad que arriben al menos 750 vehículos antes de las 72 horas. * Probabilidad Teórica => 0.08097848535892749 * Probabilidad Simulada => 0.084
A partir del generador de número al azar implementado en el ejercicio 1, y del dataset provisto, obtenido del sitio de datos abiertos del Gobierno de la Ciudad de Buenos Aires (data.buenosaires.gob.ar), el cual contiene información geográfica de barrios de la Ciudad de Buenos Aires, se pide:
class EJ5:
def __init__(self, gdf, lxm):
self.gdf = gdf.copy()
self.lxm = lxm
self.puntos_en_barrios = {}
for barrio in self.gdf['BARRIO'].tolist():
self.puntos_en_barrios[barrio] = [[],[]]
def emitir_barrios(self):
display(self.gdf['BARRIO'])
#metodo auxiliar para trasladar un barrio al primer cuadrante
def trasladar_barrio(self, datos_barrio):
extremos_barrio = datos_barrio.geometry.bounds
signo1 = 1
signo2 = 1
if(extremos_barrio.iloc[0]['minx'] < 0):
signo1 = -1
if(extremos_barrio.iloc[0]['miny'] < 0):
signo2 = -1
return datos_barrio.translate(extremos_barrio.iloc[0]['minx']*signo1, extremos_barrio.iloc[0]['miny']*signo2)
#Generar n puntos en un barrio
def generar_en_barrio(self, barrio, cantidad):
x, y = [], []
#obtengo la geometria del barrio, la traslado al primer cuadrante
#y busco el maximo entre sus extremos para formar un cuadrado que los contenga
datos_barrio = self.gdf[self.gdf['BARRIO']==barrio]
g_barrio = datos_barrio['geometry']
g_barrio_t = self.trasladar_barrio(datos_barrio)
extremos_b = datos_barrio.bounds.to_numpy()
extremos_t = g_barrio_t.bounds.to_numpy()
maximo_t = np.amax(extremos_t)
#genero puntos dentro del cuadrado que obtuve, y los acepto si están dentro del barrio
x = []
y = []
while len(x) < cantidad:
z1 = next(self.lxm) * maximo_t
z2 = next(self.lxm) * maximo_t
ps = shapely.geometry.Point(z1, z2)
if g_barrio_t.contains(ps).any():
x.append(z1 + datos_barrio.geometry.bounds.minx)
y.append(z2 + datos_barrio.geometry.bounds.miny)
self.puntos_en_barrios[barrio] = [x, y]
#Generar n puntos en todos los barrios
def generar_en_todos(self, n):
for barrio in self.gdf['BARRIO'].tolist():
self.generar_en_barrio(barrio, n)
#Graficar un barrio y los puntos generados en el mismo
def graficar_en_barrio(self, barrio):
x = self.puntos_en_barrios[barrio][0]
y = self.puntos_en_barrios[barrio][1]
title_string = "Puntos generados en el barrio de {} (n = {})".format(barrio, len(x))
fig, ax = plt.subplots(figsize=(10, 10))
fig.tight_layout(pad=3)
ax.set_title(title_string, fontsize=20)
datos_barrio = self.gdf.loc[self.gdf['BARRIO']==barrio]
datos_barrio.plot(ax=ax, color='white', edgecolor='black')
ax.scatter(x, y, s=3, color="red")
ax.grid()
return fig, ax
#Graficar el mapa completo y los puntos generados en un barrio
def graficar_mapa_completo_barrio(self, nombre_barrio):
x = self.puntos_en_barrios[nombre_barrio][0]
y = self.puntos_en_barrios[nombre_barrio][1]
title_string = "Puntos generados en el barrio de {} (n = {})".format(nombre_barrio, len(x))
fig, ax = plt.subplots(figsize=(10, 10))
fig.tight_layout(pad=3)
ax.set_title(title_string, fontsize=20)
self.gdf.plot(ax=ax, color='white', edgecolor='black')
ax.scatter(x, y, s=3, color="red")
ax.grid()
return fig, ax
#Graficar el mapa completo y los puntos generados en todos los barrios (si hay)
def graficar_mapa_completo_todos(self):
title_string = "Puntos generados en todos los barrios."
fig, ax = plt.subplots(figsize=(10, 10))
fig.tight_layout(pad=3)
ax.set_title(title_string, fontsize=20)
self.gdf.plot(ax=ax, color='white', edgecolor='black')
for b in self.gdf['BARRIO'].tolist():
x = self.puntos_en_barrios[b][0]
y = self.puntos_en_barrios[b][1]
ax.scatter(x, y, s=3, color="red")
ax.grid()
return fig, ax
def help(self):
display("generar_en_barrio('barrio', cantidad)")
display("graficar_mapa_completo_barrio('barrio')")
display("graficar_mapa_completo_todos()")
display("graficar_en_barrio('barrio')")
display("generar_en_todos(cantidad)")
ej5test = EJ5(geodataframe, LXM(42))
ej5test.generar_en_barrio('PUERTO MADERO', 100)
ej5test.generar_en_todos(100)
ej5test.graficar_en_barrio('PUERTO MADERO')
(<Figure size 1000x1000 with 1 Axes>,
<AxesSubplot: title={'center': 'Puntos generados en el barrio de PUERTO MADERO (n = 100)'}>)
ej5test.graficar_mapa_completo_barrio('PUERTO MADERO')
(<Figure size 1000x1000 with 1 Axes>,
<AxesSubplot: title={'center': 'Puntos generados en el barrio de PUERTO MADERO (n = 100)'}>)
ej5test.graficar_mapa_completo_todos()
(<Figure size 1000x1000 with 1 Axes>,
<AxesSubplot: title={'center': 'Puntos generados en todos los barrios.'}>)
# Ahora, rehagamos partes del trabajo, pero basados en un GCL en vez de en LXM
# Son los resultados distintos?
padrones = [97568, 100029, 100033, 104030]
semilla = int(np.mean(padrones))
def GCL(x, a = 1013904223, c = 1664525, m = 232):
while True:
x = ((a * x + c) % m) / m
yield x
n = 5
g = GCL(semilla)
print(f"{n} números generados al azar, en el intervalo 0-1:")
for r, i in zip(g,range(n)):
print(f"* {r}")
5 números generados al azar, en el intervalo 0-1: * 0.9224137931034483 * 0.29819411399035617 * 0.5019334719098848 * 0.895838293535956 * 0.008367325211393422
scatterplot.show()
Chi2(numero_de_muestras=10000, g=GCL(semilla)).realizar_prueba()
Frecuencias: [1005, 996, 1013, 1018, 1002, 992, 987, 1020, 994, 973] Estadístico: 1.94 Límite superior: 11.07 El test acepta la hipotesis nula
limite = 24 * 31 # Un mes
arribos_gcl = generar_poisson(limite, tasa, GCL(semilla))
print(f"Simulé {len(arribos_gcl)} arribos en {sum(arribos_gcl)} horas")
print(f"Mi tasa de arribos simulada es de {len(arribos_gcl)/sum(arribos_gcl)} arribos por hora")
comparar_poissones([(arribos_real, "Poisson Empirica"), (arribos, "Poisson Simulada con LXM"), (arribos_gcl, "Poisson Simulada con GCL")])
Simulé 7394 arribos en 744.0969199691178 horas Mi tasa de arribos simulada es de 9.936877578134409 arribos por hora
ej5test = EJ5(geodataframe, GCL(semilla))
ej5test.generar_en_barrio('PUERTO MADERO', 100)
ej5test.generar_en_todos(100)
ej5test.graficar_en_barrio('PUERTO MADERO')
(<Figure size 1000x1000 with 1 Axes>,
<AxesSubplot: title={'center': 'Puntos generados en el barrio de PUERTO MADERO (n = 100)'}>)
ej5test.graficar_mapa_completo_barrio('PUERTO MADERO')
(<Figure size 1000x1000 with 1 Axes>,
<AxesSubplot: title={'center': 'Puntos generados en el barrio de PUERTO MADERO (n = 100)'}>)
ej5test.graficar_mapa_completo_todos()
(<Figure size 1000x1000 with 1 Axes>,
<AxesSubplot: title={'center': 'Puntos generados en todos los barrios.'}>)
Conclusión: Los resultados con LXM y con el GCL a lo largo del TP nos parecen bastante similares, lo cual nos parece que tiene muchísimo sentido! LXM se basa fuertemente en un GCL.
import os
from bs4 import BeautifulSoup
wd = os.getcwd()
if 'TP2' not in wd: wd += '/TP2'
DIR = f'{wd}/paginas'
PAGINAS = os.listdir(DIR)
links = [] # Lista de tuplas de src->dst
palabras = {} # Dict de pagina: palabras
for f in PAGINAS:
with open(os.path.join(DIR, f)) as file:
soup = BeautifulSoup(file.read())
palabras[f] = soup.get_text().split()
for link in soup.find_all('a'):
dst = link.get('href').replace('http://','')
links.append((f,dst))
links
[('jonvoight.html', 'angelinajolie.html'),
('jonvoight.html', 'angelinajolie.html'),
('jonvoight.html', 'bradpitt.html'),
('bradpitt.html', 'jenniferaniston.html'),
('bradpitt.html', 'angelinajolie.html'),
('bradpitt.html', 'martinscorcese.html'),
('bradpitt.html', 'angelinajolie.html'),
('robertdeniro.html', 'martinscorcese.html'),
('angelinajolie.html', 'jonvoight.html'),
('angelinajolie.html', 'bradpitt.html')]
import networkx as nx
G = nx.DiGraph()
G.add_edges_from(links)
nx.draw_shell(G, with_labels=True, node_color='white', edge_color='grey')
matriz_ady = nx.to_pandas_adjacency(G)
matriz_ady
| jonvoight.html | angelinajolie.html | bradpitt.html | jenniferaniston.html | martinscorcese.html | robertdeniro.html | |
|---|---|---|---|---|---|---|
| jonvoight.html | 0.0 | 1.0 | 1.0 | 0.0 | 0.0 | 0.0 |
| angelinajolie.html | 1.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 |
| bradpitt.html | 0.0 | 1.0 | 0.0 | 1.0 | 1.0 | 0.0 |
| jenniferaniston.html | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| martinscorcese.html | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| robertdeniro.html | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 |
suma_filas = matriz_ady.sum(axis=1)
matriz_page_rank = matriz_ady.div(suma_filas, axis=0)
matriz_page_rank.fillna(1/len(PAGINAS), inplace=True)
matriz_page_rank
| jonvoight.html | angelinajolie.html | bradpitt.html | jenniferaniston.html | martinscorcese.html | robertdeniro.html | |
|---|---|---|---|---|---|---|
| jonvoight.html | 0.000000 | 0.500000 | 0.500000 | 0.000000 | 0.000000 | 0.000000 |
| angelinajolie.html | 0.500000 | 0.000000 | 0.500000 | 0.000000 | 0.000000 | 0.000000 |
| bradpitt.html | 0.000000 | 0.333333 | 0.000000 | 0.333333 | 0.333333 | 0.000000 |
| jenniferaniston.html | 0.166667 | 0.166667 | 0.166667 | 0.166667 | 0.166667 | 0.166667 |
| martinscorcese.html | 0.166667 | 0.166667 | 0.166667 | 0.166667 | 0.166667 | 0.166667 |
| robertdeniro.html | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 1.000000 | 0.000000 |
import numpy as np
# Hacemos N iteraciones de la matriz
iteraciones = np.linalg.matrix_power(matriz_page_rank, 50)
page_rank = {p: iteraciones[0][i] for i, p in enumerate(matriz_page_rank.index)}
page_rank
{'jonvoight.html': 0.16216216216216203,
'angelinajolie.html': 0.21621621621621548,
'bradpitt.html': 0.2432432432432427,
'jenniferaniston.html': 0.13513513513513484,
'martinscorcese.html': 0.18918918918918876,
'robertdeniro.html': 0.05405405405405395}
import random
BUSQUEDAS = 20
todas_palabras = set()
for p in palabras.values():
todas_palabras.update(p)
for p in random.choices(list(todas_palabras), k=BUSQUEDAS):
matches = [f for f in PAGINAS if p in palabras[f]]
print(f"'{p}' aparece en: {sorted(matches, key=lambda x: -page_rank[x])}")
'Twelve' aparece en: ['bradpitt.html'] 'Jolie,' aparece en: ['bradpitt.html'] '(1999).' aparece en: ['angelinajolie.html'] 'Derailed.' aparece en: ['jenniferaniston.html'] '(born' aparece en: ['bradpitt.html', 'angelinajolie.html', 'martinscorcese.html', 'jonvoight.html', 'jenniferaniston.html', 'robertdeniro.html'] 'Break-Up,' aparece en: ['jenniferaniston.html'] 'years,' aparece en: ['jonvoight.html'] 'productions' aparece en: ['bradpitt.html'] 'nomination.' aparece en: ['bradpitt.html', 'jonvoight.html'] 'beautiful' aparece en: ['angelinajolie.html'] 'named' aparece en: ['bradpitt.html'] 'four' aparece en: ['bradpitt.html'] 'Robert' aparece en: ['martinscorcese.html', 'robertdeniro.html'] 'Part' aparece en: ['robertdeniro.html'] 'Bull,' aparece en: ['robertdeniro.html'] 'Mario' aparece en: ['robertdeniro.html'] 'Leprechaun' aparece en: ['jenniferaniston.html'] 'Taxi' aparece en: ['martinscorcese.html', 'robertdeniro.html'] 'Came' aparece en: ['jenniferaniston.html'] 'stock.' aparece en: ['martinscorcese.html']
En base a lo presentado en el trabajo “Queuing theory application in imaging service analysis for small Earth observation satellites”, simular los resultados obtenidos sobre la longitud de imágenes en cola esperando ser procesadas en la sección 3.1. Pure image capture service system El ejercicio se puede resolver utilizando simpy o programación tradicional (a elección del grupo)
El paper trata de la aplicación de la teoría de colas al procesamiento de solicitudes de captura de imágenes satelitales. En particular el modelo descrito consta de dos etapas: la captura de imágenes y la descarga de imágenes. A continuación intentaremos replicar los resultados obtenidos al modelar el servicio que procesa solicitudes de captura de imágenes, descrito en la sección 3.1 del paper.
El servicio de captura de imágenes inicia con al llegada de una solicitud, y termina con la captura de la imágen solicitada.
Las solicitudes de captura de imágenes ingresan a la cola de acuerdo con un proceso de poisson de tasa $\lambda$.
Adicionalmente el modelo supone lo siguiente:
class Metricas():
def __init__(self, tasa_arribos, tasa_servicio):
self.tasa_arribos = tasa_arribos
self.tasa_servicio = tasa_servicio
self.arribos = []
self.salidas = []
self.tiempo_en_cola = []
self.tiempos_totales = []
self.cantidad_en_cola = []
self.tiempos_de_servicio = []
self.tiempos_de_servicio_cola = []
self.tasa_servicio_media_segun_largo_cola = []
self.longitudMediaTeorica = 0
self.longitudMediaEmpirica = 0
self.tiempoMedioTeorico = 0
self.tiempoMedioEmpirico = 0
def calcular_metricas_paper(self):
self.longitudMediaEmpirica = np.mean(self.cantidad_en_cola)
self.longitudMediaTeorica = (self.tasa_arribos/self.tasa_servicio)*np.e**(self.tasa_arribos/self.tasa_servicio)
self.tiempoMedioEmpirico = np.mean(self.tiempos_totales)
self.tiempoMedioTeorico = np.e**(self.tasa_arribos/self.tasa_servicio)/(self.tasa_servicio)
#print(f'Tiempo de espera promedio obtenido de la simulación: {self.tiempoMedioEmpirico}')
#print(f'Tiempo de espera promedio según el modelo: {self.tiempoMedioTeorico}')
#print(f'Longitud de la cola promedio obtenida de la simulación: {self.longitudMediaEmpirica}')
#print(f'Longitud de la cola promedio según el modelo: {self.longitudMediaTeorica}')
#generador para modelar la llegada de solicitudes
def solicitud(env, servicio, cant_arribos, tasa_arribos, tasa_servicio, metricas):
solID = 0
while solID < cant_arribos:
t_solicitud_media = 1/tasa_arribos
#t_solicitud = np.random.exponential(t_solicitud_media)
t_solicitud = generar_exponencial(tasa_arribos)
yield env.timeout(t_solicitud)
metricas.arribos.append(env.now)
solID += 1
env.process(captura_imagen(env, servicio, solID, tasa_servicio, metricas))
#generador para modelar el servicio de captura de imágenes
def captura_imagen(env, servicio, solID, tasa_servicio, metricas):
with servicio.request(priority=1) as req:
m = len(servicio.queue)
n = max(1, m)
media_tiempos_de_servicio = 1/(tasa_servicio*n)
tiempo_de_servicio = generar_exponencial(tasa_servicio)
tiempo_de_servicio_cola = tiempo_de_servicio/n
metricas.tasa_servicio_media_segun_largo_cola.append(1/tiempo_de_servicio_cola)
metricas.tiempos_de_servicio.append(tiempo_de_servicio)
metricas.tiempos_de_servicio_cola.append(tiempo_de_servicio_cola)
metricas.cantidad_en_cola.append(len(servicio.queue))
inicio_espera = env.now
yield req
fin_espera = env.now
inicio_servicio = env.now
yield env.timeout(tiempo_de_servicio_cola)
fin_servicio = env.now
metricas.tiempo_en_cola.append(fin_espera - inicio_espera)
metricas.tiempos_totales.append(fin_servicio - inicio_espera)
metricas.salidas.append(env.now)
#inicio la simulación
def simular(cant_satelites, cant_arribos, tasa_arribos, tasa_servicio, semilla):
metricas = Metricas(tasa_arribos, tasa_servicio)
env = simpy.Environment()
servicio = simpy.PriorityResource(env, capacity = cant_satelites)
np.random.seed(seed=semilla)
print(f'Simulando {cant_arribos} arribos con lambdaArribos = {tasa_arribos} y lambdaServicio = {tasa_servicio}')
env.process(solicitud(env, servicio, cant_arribos, tasa_arribos, tasa_servicio, metricas))
env.run()
return metricas
metricasHist = []
a_arribos = np.linspace(0.5, 1.5, 3)
a_llegadas = np.linspace(0.5, 1.5, 3)
longitudes = []
tiempos = []
cantidades_y_tasa_media = []
k=0
for i in a_arribos:
for j in a_llegadas:
print(k)
k=k+1
metricas = simular(cant_satelites=1, cant_arribos=100, tasa_arribos=i, tasa_servicio=j, semilla=97568)
metricas.calcular_metricas_paper()
longitudes.append((metricas.longitudMediaTeorica, metricas.longitudMediaEmpirica))
tiempos.append((metricas.tiempoMedioTeorico, metricas.tiempoMedioEmpirico))
cantidades_y_tasa_media.append((np.mean(metricas.cantidad_en_cola), (np.mean(metricas.tasa_servicio_media_segun_largo_cola))))
metricasHist.append(metricas)
0 Simulando 100 arribos con lambdaArribos = 0.5 y lambdaServicio = 0.5 1 Simulando 100 arribos con lambdaArribos = 0.5 y lambdaServicio = 1.0 2 Simulando 100 arribos con lambdaArribos = 0.5 y lambdaServicio = 1.5 3 Simulando 100 arribos con lambdaArribos = 1.0 y lambdaServicio = 0.5 4 Simulando 100 arribos con lambdaArribos = 1.0 y lambdaServicio = 1.0 5 Simulando 100 arribos con lambdaArribos = 1.0 y lambdaServicio = 1.5 6 Simulando 100 arribos con lambdaArribos = 1.5 y lambdaServicio = 0.5 7 Simulando 100 arribos con lambdaArribos = 1.5 y lambdaServicio = 1.0 8 Simulando 100 arribos con lambdaArribos = 1.5 y lambdaServicio = 1.5
fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(20, 5))
ax[0].hist(metricasHist[3].cantidad_en_cola, 20,rwidth=0.9)
ax[1].hist(metricasHist[4].cantidad_en_cola, 20,rwidth=0.9)
ax[2].hist(metricasHist[5].cantidad_en_cola, 20,rwidth=0.9)
ax[0].set_xticks(np.arange(0, 20, 1))
ax[1].set_xticks(np.arange(0, 20, 1))
ax[2].set_xticks(np.arange(0, 20, 1))
ax[0].set_xlabel(f'Cantidad de solicitudes en cola \n (λ = {metricasHist[3].tasa_arribos}, μ = {metricasHist[3].tasa_servicio})', fontsize=15)
ax[1].set_xlabel(f'Cantidad de solicitudes en cola \n (λ = {metricasHist[4].tasa_arribos}, μ = {metricasHist[4].tasa_servicio})', fontsize=15)
ax[2].set_xlabel(f'Cantidad de solicitudes en cola \n (λ = {metricasHist[5].tasa_arribos}, μ = {metricasHist[5].tasa_servicio})', fontsize=15)
fig.suptitle("Histogramas de solicitudes en cola", fontsize=20)
Text(0.5, 0.98, 'Histogramas de solicitudes en cola')
fig, ax3 = plt.subplots(nrows=1, ncols=3, figsize=(20, 5))
ax3[0].step(x=metricasHist[3].arribos, y=range(len(metricasHist[3].arribos)))
ax3[0].step(x=metricasHist[3].salidas, y=range(len(metricasHist[3].salidas)))
ax3[1].step(x=metricasHist[4].arribos, y=range(len(metricasHist[4].arribos)))
ax3[1].step(x=metricasHist[4].salidas, y=range(len(metricasHist[4].salidas)))
ax3[2].step(x=metricasHist[5].arribos, y=range(len(metricasHist[5].arribos)))
ax3[2].step(x=metricasHist[5].salidas, y=range(len(metricasHist[5].salidas)))
fig.suptitle("Evolución de arribos y salidas", fontsize=20)
ax3[0].set_xlabel(f'Arribos y salidas \n (λ = {metricasHist[3].tasa_arribos}, μ = {metricasHist[3].tasa_servicio})', fontsize=15)
ax3[1].set_xlabel(f'Arribos y salidas \n (λ = {metricasHist[4].tasa_arribos}, μ = {metricasHist[4].tasa_servicio})', fontsize=15)
ax3[2].set_xlabel(f'Arribos y salidas \n (λ = {metricasHist[5].tasa_arribos}, μ = {metricasHist[5].tasa_servicio})', fontsize=15)
ax3[0].legend([f'λ = {metricasHist[3].tasa_arribos}, μ = {metricasHist[3].tasa_servicio}'])
ax3[1].legend([f'λ = {metricasHist[4].tasa_arribos}, μ = {metricasHist[4].tasa_servicio}'])
ax3[2].legend([f'λ = {metricasHist[5].tasa_arribos}, μ = {metricasHist[5].tasa_servicio}'])
<matplotlib.legend.Legend at 0x7f77a89bb8e0>
fig, ax4 = plt.subplots(nrows=2, ncols=1, figsize=(20, 10))
p = sns.kdeplot(metricasHist[0].cantidad_en_cola,ax=ax4[0])
p = sns.kdeplot(metricasHist[3].cantidad_en_cola,ax=ax4[0])
p = sns.kdeplot(metricasHist[6].cantidad_en_cola,ax=ax4[0])
#p = sns.kdeplot(metricasHist[12].cantidad_en_cola,ax=ax4[0])
p = sns.kdeplot(metricasHist[0].cantidad_en_cola,ax=ax4[1])
p = sns.kdeplot(metricasHist[1].cantidad_en_cola,ax=ax4[1])
p = sns.kdeplot(metricasHist[2].cantidad_en_cola,ax=ax4[1])
#p = sns.kdeplot(metricasHist[3].cantidad_en_cola,ax=ax4[1])
ax4[0].set_xticks(np.arange(0, 30, 1))
ax4[0].legend([f'λ = {metricasHist[0].tasa_arribos}, μ = {metricasHist[0].tasa_servicio}',
f'λ = {metricasHist[3].tasa_arribos}, μ = {metricasHist[3].tasa_servicio}',
f'λ = {metricasHist[6].tasa_arribos}, μ = {metricasHist[6].tasa_servicio}'])
ax4[0].set_xlabel('Solicitudes en cola', fontsize=15)
ax4[0].set_ylabel('Densidad', fontsize=15)
ax4[0].set_title("Densidad de cantidad de solicitudes en cola", fontsize=20)
ax4[1].set_xticks(np.arange(0, 30, 1))
ax4[1].legend([f'λ = {metricasHist[0].tasa_arribos}, μ = {metricasHist[0].tasa_servicio}',
f'λ = {metricasHist[1].tasa_arribos}, μ = {metricasHist[1].tasa_servicio}',
f'λ = {metricasHist[2].tasa_arribos}, μ = {metricasHist[2].tasa_servicio}'])
ax4[1].set_xlabel('Solicitudes en cola', fontsize=15)
ax4[1].set_ylabel('Densidad', fontsize=15)
ax4[1].set_title("Densidad de cantidad de solicitudes en cola", fontsize=20)
fig.tight_layout()
fig, ax5 = plt.subplots(nrows=1, ncols=1, figsize=(20, 5))
ax5.scatter(*zip(*cantidades_y_tasa_media))
ax5.set_xlabel('Mean image capture queue length', fontsize=15)
ax5.set_ylabel('Mean image capture service rate (req/day)', fontsize=15)
ax5.set_title("Longitudes de cola medias", fontsize=20)
Text(0.5, 1.0, 'Longitudes de cola medias')
fig, ax6 = plt.subplots(nrows=1, ncols=2, figsize=(20, 5))
ax6[0].scatter(*zip(*longitudes))
ax6[0].set_xlabel('Longitud Media Teórica', fontsize=15)
ax6[0].set_ylabel('Longitud Media Empirica', fontsize=15)
ax6[0].set_title("Longitudes de cola medias", fontsize=20)
ax6[1].scatter(*zip(*tiempos))
ax6[1].set_xlabel('Tiempo medio de servicio Teórico', fontsize=15)
ax6[1].set_ylabel('Tiempo medio de servicio Empírico', fontsize=15)
ax6[1].set_title("Tiempos medios de servicio", fontsize=20)
Text(0.5, 1.0, 'Tiempos medios de servicio')
fig, ax7 = plt.subplots(nrows=3, ncols=2, figsize=(20,10))
p = sns.kdeplot(metricasHist[3].tiempos_de_servicio,ax=ax7[0,0])
p = sns.kdeplot(metricasHist[3].tiempos_de_servicio_cola,ax=ax7[0,1])
p = sns.kdeplot(metricasHist[4].tiempos_de_servicio,ax=ax7[1,0])
p = sns.kdeplot(metricasHist[4].tiempos_de_servicio_cola,ax=ax7[1,1])
p = sns.kdeplot(metricasHist[5].tiempos_de_servicio,ax=ax7[2,0])
p = sns.kdeplot(metricasHist[5].tiempos_de_servicio_cola,ax=ax7[2,1])
ax7[0,0].legend([f'λ = {metricasHist[3].tasa_arribos}, μ = {metricasHist[3].tasa_servicio}'])
ax7[0,1].legend([f'λ = {metricasHist[3].tasa_arribos}, μ = {metricasHist[3].tasa_servicio}'])
ax7[1,0].legend([f'λ = {metricasHist[4].tasa_arribos}, μ = {metricasHist[4].tasa_servicio}'])
ax7[1,1].legend([f'λ = {metricasHist[4].tasa_arribos}, μ = {metricasHist[4].tasa_servicio}'])
ax7[2,0].legend([f'λ = {metricasHist[5].tasa_arribos}, μ = {metricasHist[5].tasa_servicio}'])
ax7[2,1].legend([f'λ = {metricasHist[5].tasa_arribos}, μ = {metricasHist[5].tasa_servicio}'])
ax7[0,0].set_yticks(np.arange(0, 1.6, 0.2))
ax7[0,1].set_yticks(np.arange(0, 1.6, 0.2))
ax7[1,0].set_yticks(np.arange(0, 1.6, 0.2))
ax7[1,1].set_yticks(np.arange(0, 1.6, 0.2))
ax7[2,0].set_yticks(np.arange(0, 1.6, 0.2))
ax7[2,1].set_yticks(np.arange(0, 1.6, 0.2))
ax7[0,1].legend([f'λ = {metricasHist[3].tasa_arribos}, μ = {metricasHist[3].tasa_servicio}'])
ax7[1,0].legend([f'λ = {metricasHist[4].tasa_arribos}, μ = {metricasHist[4].tasa_servicio}'])
ax7[1,1].legend([f'λ = {metricasHist[4].tasa_arribos}, μ = {metricasHist[4].tasa_servicio}'])
ax7[2,0].legend([f'λ = {metricasHist[5].tasa_arribos}, μ = {metricasHist[5].tasa_servicio}'])
ax7[2,1].legend([f'λ = {metricasHist[5].tasa_arribos}, μ = {metricasHist[5].tasa_servicio}'])
#fig.tight_layout()
fig.suptitle("Densidad de tiempos de servicio", fontsize=20)
Text(0.5, 0.98, 'Densidad de tiempos de servicio')
from scipy.stats import expon, uniform
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
from math import floor
sns.set(style="darkgrid")
params = {
'axes.titlesize': 20,
'axes.labelsize': 15,
'legend.fontsize': 20,
'figure.figsize': (10, 6),
}
pylab.rcParams.update(params)
class Consulta:
def __init__(self, tiempo_arribo, duracion, tiempo_finalizacion_anterior):
if tiempo_arribo < tiempo_finalizacion_anterior:
self.duracion_espera = tiempo_finalizacion_anterior - tiempo_arribo
self.tuvo_que_esperar = True
else:
self.duracion_espera = 0
self.tuvo_que_esperar = False
self.tiempo_finalizacion = tiempo_arribo + duracion + self.duracion_espera
def simular_consultas(tiempos_de_arribo, duraciones):
consultas = []
# Uso -1 como tiempo de finalizacion del anterior para que empieze ni bien llega
primera_consulta = Consulta(tiempos_de_arribo[0], duraciones[0], tiempo_finalizacion_anterior = -1)
consultas.append(primera_consulta)
for tiempo_arribo, duracion in zip(tiempos_de_arribo[1:], duraciones[1:]):
consulta = Consulta(tiempo_arribo,
duracion,
tiempo_finalizacion_anterior = consultas[-1].tiempo_finalizacion)
consultas.append(consulta)
tiempo_espera_total = 0
cant_consultas_que_no_esperaron = 0
for consulta in consultas:
tiempo_espera_total += consulta.duracion_espera
if not consulta.tuvo_que_esperar:
cant_consultas_que_no_esperaron += 1
return tiempo_espera_total, cant_consultas_que_no_esperaron
def mostrar_resultados(titulo, tiempo_espera_total, cant_consultas_que_no_esperaron, cant_consultas):
print(f"--------------{titulo}--------------")
segundos = round(tiempo_espera_total / cant_consultas, 0)
minutos = floor(segundos/60)
if(minutos != 0):
segundos = segundos % (minutos * 60)
horas = floor(minutos/60)
if(horas != 0):
minutos = minutos % (horas * 60)
dias = floor(horas/24)
if(dias != 0):
horas = horas % (dias * 24)
print(f"\tEl tiempo de espera promedio fue ", end="")
if(dias != 0):
print(str(dias) + " Dias ", end="")
if(horas != 0):
print(str(horas) + " Horas ", end="")
if(minutos != 0):
print(str(minutos) + " Minutos ", end="")
if(segundos != 0):
print(str(segundos) + " Segundos")
print(f"\tEl {round(cant_consultas_que_no_esperaron/cant_consultas, 5)}% de las consultas no tuvo que esperar")
def simular_consultas_centralizado(cant_consultas, media_tiempos_de_arribo, media_duraciones_consultas):
tiempos_de_arribo = expon.rvs(size=cant_consultas, scale=media_tiempos_de_arribo, random_state=10)
duraciones = expon.rvs(size=cant_consultas, scale=media_duraciones_consultas, random_state=10)
return simular_consultas(tiempos_de_arribo, duraciones)
t, f = simular_consultas_centralizado(100000, 4, 0.8)
mostrar_resultados("Base de datos central", t,f, 100000)
--------------Base de datos central-------------- El tiempo de espera promedio fue 11 Horas 3 Minutos 16.0 Segundos El 2e-05% de las consultas no tuvo que esperar
def simular_consultas_distribuido(cant_consultas, media_tiempos_de_arribo,
media_duraciones_consultas1, media_duraciones_consultas2, p):
tiempos_de_arribo = expon.rvs(size=cant_consultas, scale=media_tiempos_de_arribo, random_state=10)
u = uniform.rvs(size=cant_consultas, random_state=10)
tiempos_de_arribo1 = []
tiempos_de_arribo2 = []
for ui, tiempo_de_arribo in zip(u, tiempos_de_arribo):
if p <= ui:
tiempos_de_arribo1.append(tiempo_de_arribo)
else:
tiempos_de_arribo2.append(tiempo_de_arribo)
duraciones1 = expon.rvs(size=len(tiempos_de_arribo1), scale=media_duraciones_consultas1, random_state=10)
t1, f1 = simular_consultas(tiempos_de_arribo1, duraciones1)
duraciones2 = expon.rvs(size=len(tiempos_de_arribo2), scale=media_duraciones_consultas2, random_state=10)
t2, f2 = simular_consultas(tiempos_de_arribo2, duraciones2)
return t1, f1, len(tiempos_de_arribo1), t2, f2, len(tiempos_de_arribo2)
t1, f1, l1, t2, f2, l2 = simular_consultas_distribuido(100000, 4, 0.7, 1, 0.7)
mostrar_resultados("Ambas bases distribuidas", t1 + t2, f2 + f2, l1 + l2)
mostrar_resultados("Base de datos distribuida 1", t1, f1, l1)
mostrar_resultados("Base de datos distribuida 2", t2, f2, l2)
--------------Ambas bases distribuidas-------------- El tiempo de espera promedio fue 7 Horas 40 Minutos 21.0 Segundos El 4e-05% de las consultas no tuvo que esperar --------------Base de datos distribuida 1-------------- El tiempo de espera promedio fue 2 Horas 51 Minutos 34.0 Segundos El 7e-05% de las consultas no tuvo que esperar --------------Base de datos distribuida 2-------------- El tiempo de espera promedio fue 9 Horas 42 Minutos 30.0 Segundos El 3e-05% de las consultas no tuvo que esperar
def simular_variaciones(cantidades_consultas=[100000],
medias_tiempos_de_arribo=[4],
medias_duraciones_consultas = [0.8],
medias_duraciones_consultas1 = [0.7],
medias_duraciones_consultas2 = [1],
lista_p = [0.7]):
tiempos_centralizado = []
frecuencias_centralizado = []
tiempos_distribuido = []
frecuencias_distribuido = []
for cant_consultas in cantidades_consultas:
for media_tiempos_de_arribo in medias_tiempos_de_arribo:
for media_duraciones_consultas in medias_duraciones_consultas:
for media_duraciones_consultas1 in medias_duraciones_consultas1:
for media_duraciones_consultas2 in medias_duraciones_consultas2:
for p in lista_p:
t, f = simular_consultas_centralizado(cant_consultas, media_tiempos_de_arribo,
media_duraciones_consultas)
tiempos_centralizado.append(t / cant_consultas)
frecuencias_centralizado.append(f /cant_consultas)
t1, f1, l1, t2, f2, l2 = simular_consultas_distribuido(cant_consultas,
media_tiempos_de_arribo,
media_duraciones_consultas1,
media_duraciones_consultas2,
p)
tiempos_distribuido.append((t1 + t2) / (l1 + l2))
frecuencias_distribuido.append((f1 + f2) / (l1 + l2))
return tiempos_centralizado, frecuencias_centralizado, tiempos_distribuido, frecuencias_distribuido
def graficar_tiempo_espera(x, xlabel, tiempos_centralizado, tiempos_distribuido,
centralizado_es_referencia = False, distribuido_es_referencia = False):
if centralizado_es_referencia:
color_centralizado = 'b'
linestyle_centralizado = '--'
else:
color_centralizado = 'bo'
linestyle_centralizado = '-'
if distribuido_es_referencia:
color_distribuido = 'g'
linestyle_distribuido = '--'
else:
color_distribuido = 'go'
linestyle_distribuido = '-'
plt.plot(x, tiempos_centralizado, color_centralizado, linestyle=linestyle_centralizado, label="Base centralizada")
plt.plot(x, tiempos_distribuido, color_distribuido, linestyle=linestyle_distribuido, label="Bases distribuidas")
plt.title(f"Tiempo de espera medio según\n{xlabel.lower()}")
plt.xlabel(xlabel)
plt.ylabel("Tiempo de espera medio")
plt.legend()
plt.show()
def graficar_fraccion_que_no_espero(x, xlabel, tiempos_centralizado, tiempos_distribuido,
centralizado_es_referencia = False, distribuido_es_referencia = False):
if centralizado_es_referencia:
color_centralizado = 'b'
linestyle_centralizado = '--'
else:
color_centralizado = 'bo'
linestyle_centralizado = '-'
if distribuido_es_referencia:
color_distribuido = 'g'
linestyle_distribuido = '--'
else:
color_distribuido = 'go'
linestyle_distribuido = '-'
plt.plot(x, tiempos_centralizado, color_centralizado, linestyle=linestyle_centralizado, label="Base centralizada")
plt.plot(x, tiempos_distribuido, color_distribuido, linestyle=linestyle_distribuido, label="Bases distribuidas")
plt.title(f"Fracción de solicitudes que no esperaron según\n{xlabel.lower()}")
plt.xlabel(xlabel)
plt.ylabel("Fracción de solicitudes que no esperaron")
plt.legend()
medias_tiempos_de_arribo = range(1, 11)
tiempos_centralizado, frecuencias_centralizado, tiempos_distribuido, frecuencias_distribuido = simular_variaciones(
medias_tiempos_de_arribo = medias_tiempos_de_arribo
)
graficar_tiempo_espera(medias_tiempos_de_arribo, "Media de los tiempos de arribo de consultas",
tiempos_centralizado, tiempos_distribuido)
graficar_fraccion_que_no_espero(medias_tiempos_de_arribo, "Media de los tiempos de arribo de consultas",
frecuencias_centralizado, frecuencias_distribuido)
lambda_l = []
for i in range(2, 20, 2):
lambda_l.append(i/10)
tiempos_centralizado, frecuencias_centralizado, tiempos_distribuido, frecuencias_distribuido = simular_variaciones(
medias_duraciones_consultas = lambda_l
)
graficar_tiempo_espera(lambda_l, "Lambda",
tiempos_centralizado, tiempos_distribuido, distribuido_es_referencia = True)
graficar_fraccion_que_no_espero(lambda_l, "Lambda",
frecuencias_centralizado, frecuencias_distribuido, distribuido_es_referencia = True)
lambda_1 = []
for i in range(2, 20, 2):
lambda_1.append(i/10)
tiempos_centralizado, frecuencias_centralizado, tiempos_distribuido, frecuencias_distribuido = simular_variaciones(
medias_duraciones_consultas1 = lambda_1
)
graficar_tiempo_espera(lambda_1, "Lambda1",
tiempos_centralizado, tiempos_distribuido, centralizado_es_referencia = True)
graficar_fraccion_que_no_espero(lambda_1, "Lambda1",
frecuencias_centralizado, frecuencias_distribuido, centralizado_es_referencia = True)
lambda_2 = []
for i in range(2, 20, 2):
lambda_2.append(i/10)
tiempos_centralizado, frecuencias_centralizado, tiempos_distribuido, frecuencias_distribuido = simular_variaciones(
medias_duraciones_consultas1 = lambda_2
)
graficar_tiempo_espera(lambda_2, "Lambda2",
tiempos_centralizado, tiempos_distribuido, centralizado_es_referencia = True)
graficar_fraccion_que_no_espero(lambda_2, "Lambda2",
frecuencias_centralizado, frecuencias_distribuido, centralizado_es_referencia = True)
lista_p = []
for i in range(1, 10):
lista_p.append(i/10)
tiempos_centralizado, frecuencias_centralizado, tiempos_distribuido, frecuencias_distribuido = simular_variaciones(
lista_p = lista_p
)
graficar_tiempo_espera(lista_p, "p",
tiempos_centralizado, tiempos_distribuido, centralizado_es_referencia = True)
graficar_fraccion_que_no_espero(lista_p, "p",
frecuencias_centralizado, frecuencias_distribuido, centralizado_es_referencia = True)
a. Simular 1000 días completos de 24 hrs.
b. Para un día en particular graficar la cantidad de billetes en el cajero luego de cada transacción.
c. Calcular el tiempo medio que los clientes demoraron en el sistema (espera + utilización del cajero)
d. ¿Recomienda a la entidad que implemente el cambio de cajero?
from scipy.stats import expon, uniform
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
from math import floor
from random import randint
import random
import matplotlib.pyplot as plt
from statistics import median
class Cajero:
def __init__(self):
self._billetes_en_cajero = 2000 # son todos de 100, arranca el día lleno
self._capacidad_maxima = 2000 #billetes max q entran
self._historial_dia_actual = []
self._historial = []
def nuevo_dia(self):
self._billetes_en_cajero = 2000
self._historial.append(self._historial_dia_actual)
self._historial_dia_actual = []
def retirar(self, cantidad_de_billetes):
if self._billetes_en_cajero < cantidad_de_billetes:
return False
self._billetes_en_cajero -= cantidad_de_billetes
self._historial_dia_actual.append(self._billetes_en_cajero)
return True
def depositar(self, cantidad_de_billetes):
if self._billetes_en_cajero + cantidad_de_billetes > self._capacidad_maxima:
return False
self._billetes_en_cajero += cantidad_de_billetes
self._historial_dia_actual.append(self._billetes_en_cajero)
return True
def historial(self):
return self._historial
class Banco:
def __init__(self):
self._cajero = Cajero()
self._minutos = 0 # al pasa 24 * 60 termina el día y restarteamos
self._tiempos_de_espera = []
def _simular_dia(self):
print(expon.rvs(size=3, random_state=10))
# tarda en llegar n
def _grupo_cliente(self):
return random.choices([1,2], weights=(75, 25))[0]
def historial_de_cajero(self):
return self._cajero.historial()
def simular(self, dias, imprimir_estadisticas_clientes=False, imprimir_media_espera=False):
tiempos_de_arribo = expon.rvs(size=144 * 60 * dias, scale=10, random_state=10)
duraciones_1 = expon.rvs(size=len(tiempos_de_arribo), scale=1.5, random_state=10)
duraciones_2 = expon.rvs(size=len(tiempos_de_arribo), scale=5, random_state=10)
clientes_perdidos = 0
cantidad_de_clientes = 0
cantidad_de_dias = 0
cantidad_de_clientes = 0
minutos_en_hora = 24 * 60
for i, tiempo in enumerate(tiempos_de_arribo):
grupo_cliente = self._grupo_cliente()
self._minutos += tiempo
duracion_de_persona = 0
if self._minutos >= minutos_en_hora:
self._minutos = minutos_en_hora - self._minutos
self._cajero.nuevo_dia()
cantidad_de_dias += 1
if cantidad_de_dias == dias:
break
if grupo_cliente == 1: # Retiran
cant_billetes = randint(3, 50)
fue_exitoso = self._cajero.retirar(cant_billetes)
clientes_perdidos += 0 if fue_exitoso else 1
duracion_de_persona = duraciones_1[i]
else: # Depositan
cant_billetes = randint(10, 110)
fue_exitoso = self._cajero.depositar(cant_billetes)
clientes_perdidos += 0 if fue_exitoso else 1
duracion_de_persona = duraciones_2[i]
demora = 0
if i - 1 < len(tiempos_de_arribo) and tiempos_de_arribo[i+1] < duracion_de_persona:
demora = duracion_de_persona - tiempos_de_arribo[i+1]
self._minutos += demora
self._tiempos_de_espera.append(duracion_de_persona + demora)
cantidad_de_clientes += 1
if imprimir_estadisticas_clientes:
print(f'Pasaron {cantidad_de_dias} días')
print(f'Cantidad de clientes: {cantidad_de_clientes}')
print(f'Se perdieron {clientes_perdidos} clientes')
porcentaje_clientes_perdidos = clientes_perdidos * 100 / cantidad_de_clientes
print(f'Porcentaje de clientes perdidos: {floor(porcentaje_clientes_perdidos)}%')
if imprimir_media_espera:
mediana = round(median(self._tiempos_de_espera), 2)
print(f'Tiempo medio de espera: {mediana} minutos')
# Simular 1000 días
banco = Banco()
banco.simular(dias=1000, imprimir_estadisticas_clientes=True)
Pasaron 1000 días Cantidad de clientes: 138501 Se perdieron 2587 clientes Porcentaje de clientes perdidos: 1%
# cantidad de billetes en el cajero dps de cada transacción
banco = Banco()
banco.simular(dias=1000)
monto_en_cajero = banco.historial_de_cajero()[1]
# make up some data
x = [i for i,monto in enumerate(monto_en_cajero)]
# plot
plt.scatter(x,monto_en_cajero, linewidths = 0.1)
# beautify the x-labels
plt.gcf().autofmt_xdate()
# Calcular el tiempo medio que los clientes demoraron en el sistema (espera + utilización del cajero)
banco = Banco()
banco.simular(dias=1000, imprimir_media_espera=True)
Tiempo medio de espera: 1.38 minutos
Sí, ya que el promedio de clientes perdidos es entre 1% y 2%, mientras que nos indicaron antes que la perdida de clientes era de una 20%